This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.
Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!
SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.
Morrisey Method
The Morrisey Method is taken from Ma, Fassnacht and Kampf..
In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37
In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler
San Juan Area SNOTEL sites:
Beartown 327 Original 4/18/2005
Cascade #2 389 Morrisey 6/18/2004
Cumbres Trestle 431 Oyler -> Morrisey 6/8/2005
Idarado 538 Morrisey 7/12/2004
Lone Cone 589 Oyler -> Morrisey 6/22/2005
Middle Creek 624 Morrisey 6/26/2006
Mineral Creek 629 Oyler ?? -> Morrisey 6/22/2004
Molas Lake 632 NSCE-Morrisey, Bias-Original 10/2/2003
Red Mountain Pass 713 Morrisey 8/18/2004
Scotch Creek 739 Morrisey 6/15/2004
Slumgullion 762 Morrisey 6/26/2006
Spud Mountain 780 Morrisey 6/28/2004
Stump Lakes 797 Morrisey 7/22/2005
Upper Rio Grande 839 Oyler -> Morrisey 5/26/2004 *Why is this in red?
Upper San Juan 840 Morrisey 4/7/2004
Vallecito 843 Morrisey 7/22/2005
Wolf Creek Summit 874 Morrisey 7/12/2004
SNOTEL_san_juan_Area <- snotel_download(site_id = c(327, 389, 431, 538, 589, 624, 629, 632, 713, 739, 762, 780, 797, 839, 840, 843, 874), path = tempdir('../data'), internal = TRUE)
write.csv(SNOTEL_san_juan_Area,"C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", row.names = FALSE) #write in the raw data
Original
SNOTEL_san_juan_Area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", header = TRUE)
snotel_327 <- SNOTEL_san_juan_Area %>%
filter(site_id == "327")
#str(snotel_327) # check the date, usually a character.
snotel_327$Date <- as.Date(snotel_327$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_327_clean <- snotel_327 %>% # filter for the timeframe
filter(Date >= "1983-08-10" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_327_clean <- snotel_327_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_327_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_327_clean <- snotel_327_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40)
ggplot(snotel_327_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_327_cull_count <- snotel_327_clean %>%
filter(temperature_min < -40) %>%
count(waterYear)
snotel_327_cull_count
## # A tibble: 0 x 2
## # Groups: waterYear [0]
## # ... with 2 variables: waterYear <dbl>, n <int>
# filtering for too few observations in a year
snotel_327_cull_count_days <- snotel_327_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_327_cull_count_days
## # A tibble: 3 x 2
## # Groups: waterYear [3]
## waterYear n
## <dbl> <int>
## 1 1983 52
## 2 1994 337
## 3 2003 346
1983, 1994, 2003 need to be culled.
snotel_327_clean_culled <- snotel_327_clean %>%
filter(waterYear != "1983" & waterYear != "1994" & waterYear != "2003") #%>%
#filter(temperature_mean > -49)
2005-04-18 installed ext range sensor. Original
snotel_327_adjusted <- snotel_327_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-04-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_327 <- snotel_327_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_327 <- yearly_wy_aver_327 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_327 <- daily_wy_aver_327 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_327$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_327 <-daily_wy_aver_327 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_327$date_temp <- signif(daily_wy_aver2_327$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_327, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_327 <- daily_wy_aver_327 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_327 <- standard_dev_327 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_327 <- standard_dev_all_327 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_327 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1984 | 8.215357 |
| 1985 | 8.122497 |
| 1986 | 7.129599 |
| 1987 | 7.619777 |
| 1988 | 8.257318 |
| 1989 | 8.044976 |
| 1990 | 7.675963 |
| 1991 | 7.684970 |
| 1992 | 7.079517 |
| 1993 | 7.642300 |
| 1995 | 7.304999 |
| 1996 | 7.681058 |
| 1997 | 7.711950 |
| 1998 | 8.105700 |
| 1999 | 6.825322 |
| 2000 | 7.526413 |
| 2001 | 7.912097 |
| 2002 | 8.099855 |
| 2004 | 7.622635 |
| 2005 | 7.225545 |
| 2006 | 7.107924 |
| 2007 | 7.426444 |
| 2008 | 7.802192 |
| 2009 | 6.893765 |
| 2010 | 7.820347 |
| 2011 | 7.632307 |
| 2012 | 7.466342 |
| 2013 | 8.000100 |
| 2014 | 7.260505 |
| 2015 | 6.671562 |
| 2016 | 7.474460 |
| 2017 | 7.072348 |
| 2018 | 6.933787 |
| 2019 | 7.658054 |
| 2020 | 7.769482 |
| 2021 | 7.841403 |
| 2022 | 7.318492 |
ggplot(standard_dev_all_327, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 327 average temperatures for water years 2005-2021
sd_mk_327 <- mk.test(standard_dev_all_327$sd_2)
print(sd_mk_327)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -140.0000000 5846.0000000 -0.2102102
sd_sens_327 <- sens.slope(standard_dev_all_327$sd_2)
print(sd_sens_327)
##
## Sen's slope
##
## data: standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.025461316 0.001275883
## sample estimates:
## Sen's slope
## -0.01223021
Morrisey 6/18/2004
snotel_389 <- SNOTEL_san_juan_Area %>%
filter(site_id == "389")
#str(snotel_389) # check the date, usually a character.
snotel_389$Date <- as.Date(snotel_389$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_389_clean <- snotel_389 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_389_clean <- snotel_389_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_389_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_389_clean <- snotel_389_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40) %>%
filter(temperature_mean > -40)
ggplot(snotel_389_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
temp_389_xts <- xts(snotel_389_clean$temperature_mean, order.by = snotel_389_clean$Date)
dygraph(temp_389_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_389_cull_count <- snotel_389_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_389_cull_count
## # A tibble: 40 x 2
## # Groups: waterYear [40]
## waterYear n
## <dbl> <int>
## 1 1983 325
## 2 1984 350
## 3 1985 350
## 4 1986 365
## 5 1987 358
## 6 1988 366
## 7 1989 330
## 8 1990 348
## 9 1991 364
## 10 1992 366
## # ... with 30 more rows
# filtering for too few observations in a year
snotel_389_cull_count_days <- snotel_389_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_389_cull_count_days
## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1983 325
## 2 1989 330
## 3 1990 348
## 4 1994 310
## 5 2006 330
## 6 2010 334
## 7 2014 314
## 8 2015 185
## 9 2016 334
1983, 1989, 1990, 1994, 2006, 2010, 2014, 2015, 2016 need to be culled.
snotel_389_clean_culled <- snotel_389_clean %>%
filter(waterYear != "1983" & waterYear != "1990" & waterYear != "1994" & waterYear != "2006" & waterYear != "2010" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016") #%>%
#filter(temperature_mean > -49)
Morrisey 6/18/2004
snotel_389_adjusted <- snotel_389_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-06-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_389 <- snotel_389_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_389 <- yearly_wy_aver_389 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_389 <- daily_wy_aver_389 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_389$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_389 <-daily_wy_aver_389 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_389$date_temp <- signif(daily_wy_aver2_389$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_389, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_389 <- daily_wy_aver_389 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389 <- standard_dev_389 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_389 <- standard_dev_all_389 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_389 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1984 | 10.776321 |
| 1985 | 7.879611 |
| 1986 | 10.887499 |
| 1987 | 9.408888 |
| 1988 | 10.711609 |
| 1989 | 10.418326 |
| 1991 | 10.101153 |
| 1992 | 8.637142 |
| 1993 | 8.884440 |
| 1995 | 9.371322 |
| 1996 | 9.429341 |
| 1997 | 10.027003 |
| 1998 | 9.698491 |
| 1999 | 9.029327 |
| 2000 | 9.288091 |
| 2001 | 10.630687 |
| 2002 | 10.346726 |
| 2003 | 9.727202 |
| 2004 | 9.282545 |
| 2005 | 8.472227 |
| 2007 | 9.686674 |
| 2008 | 9.541070 |
| 2009 | 8.659960 |
| 2011 | 9.334867 |
| 2012 | 9.418612 |
| 2013 | 9.751020 |
| 2017 | 9.144382 |
| 2018 | 8.911243 |
| 2019 | 9.427120 |
| 2020 | 9.880776 |
| 2021 | 9.863539 |
| 2022 | 9.441803 |
ggplot(standard_dev_all_389, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 389 average temperatures for water years 2005-2021
sd_mk_389 <- mk.test(standard_dev_all_389$sd_2)
print(sd_mk_389)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -60.0000000 3802.6666667 -0.1209677
sd_sens_389 <- sens.slope(standard_dev_all_389$sd_2)
print(sd_sens_389)
##
## Sen's slope
##
## data: standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.04702986 0.01336703
## sample estimates:
## Sen's slope
## -0.01646233
#using the clean culled df:
#average water year temperature
yearly_wy_aver_389_ad <- snotel_389_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_389_ad <- yearly_wy_aver_389_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_389_ad <- daily_wy_aver_389_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_389_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_389_ad <-daily_wy_aver_389_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_389_ad$date_temp_ad <- signif(daily_wy_aver2_389_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_389_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_389_ad <- daily_wy_aver_389_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389_ad <- standard_dev_389_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_389_ad <- standard_dev_all_389_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_389_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1984 | 10.061639 |
| 1985 | 7.175032 |
| 1986 | 9.930924 |
| 1987 | 8.682951 |
| 1988 | 9.868631 |
| 1989 | 9.733860 |
| 1991 | 9.472969 |
| 1992 | 8.002025 |
| 1993 | 8.283300 |
| 1995 | 8.657379 |
| 1996 | 8.697682 |
| 1997 | 9.311693 |
| 1998 | 8.929771 |
| 1999 | 8.328247 |
| 2000 | 8.510889 |
| 2001 | 9.808164 |
| 2002 | 9.563945 |
| 2003 | 8.940614 |
| 2004 | 8.587303 |
| 2005 | 8.472294 |
| 2007 | 9.686174 |
| 2008 | 9.541322 |
| 2009 | 8.659292 |
| 2011 | 9.334237 |
| 2012 | 9.418205 |
| 2013 | 9.750812 |
| 2017 | 9.144257 |
| 2018 | 8.911176 |
| 2019 | 9.426664 |
| 2020 | 9.880403 |
| 2021 | 9.863015 |
| 2022 | 9.441185 |
ggplot(standard_dev_all_389_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 389 average temperatures for water years 1986-2021
sd_mk_389_ad <- mk.test(standard_dev_all_389_ad$sd_2)
print(sd_mk_389_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 64.0000000 3802.6666667 0.1290323
sd_sens_389_ad <- sens.slope(standard_dev_all_389_ad$sd_2)
print(sd_sens_389_ad)
##
## Sen's slope
##
## data: standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01140359 0.04896578
## sample estimates:
## Sen's slope
## 0.01785583
Oyler -> Morrisey 6/8/2005
snotel_431 <- SNOTEL_san_juan_Area %>%
filter(site_id == "431")
#str(snotel_431) # check the date, usually a character.
snotel_431$Date <- as.Date(snotel_431$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_431_clean <- snotel_431 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_431_clean <- snotel_431_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_431_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_431_clean <- snotel_431_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40) %>%
filter(temperature_min < 25)
ggplot(snotel_431_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_431_cull_count <- snotel_431_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_431_cull_count
## # A tibble: 36 x 2
## # Groups: waterYear [36]
## waterYear n
## <dbl> <int>
## 1 1987 1
## 2 1988 341
## 3 1989 364
## 4 1990 348
## 5 1991 341
## 6 1992 312
## 7 1993 356
## 8 1994 317
## 9 1995 364
## 10 1996 364
## # ... with 26 more rows
# filtering for too few observations in a year
snotel_431_cull_count_days <- snotel_431_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_431_cull_count_days
## # A tibble: 7 x 2
## # Groups: waterYear [7]
## waterYear n
## <dbl> <int>
## 1 1987 1
## 2 1988 341
## 3 1990 348
## 4 1991 341
## 5 1992 312
## 6 1994 317
## 7 2014 288
1987, 1988, 1990, 1991, 1992, 1994, 2014 need to be culled.
snotel_431_clean_culled <- snotel_431_clean %>%
filter(waterYear != "1987" & waterYear != "1988" & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_431_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)
dygraph(temp_431_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_431_clean_culled <- snotel_431_clean_culled %>%
filter(temperature_mean < 20)
temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)
dygraph(temp_431_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
Oyler -> Morrisey 6/8/2005
snotel_431_adjusted <- snotel_431_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-08", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_431 <- snotel_431_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_431 <- yearly_wy_aver_431 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_431 <- daily_wy_aver_431 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_431$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_431 <-daily_wy_aver_431 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_431$date_temp <- signif(daily_wy_aver2_431$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_431, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_431 <- daily_wy_aver_431 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431 <- standard_dev_431 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_431 <- standard_dev_all_431 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_431 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1989 | 8.484011 |
| 1993 | 8.101376 |
| 1995 | 7.750446 |
| 1996 | 7.952021 |
| 1997 | 8.217323 |
| 1998 | 8.600343 |
| 1999 | 7.226028 |
| 2000 | 8.164156 |
| 2001 | 8.669911 |
| 2002 | 9.005161 |
| 2003 | 8.459072 |
| 2004 | 8.306872 |
| 2005 | 8.304679 |
| 2006 | 7.417238 |
| 2007 | 7.937180 |
| 2008 | 8.259566 |
| 2009 | 7.296965 |
| 2010 | 8.418225 |
| 2011 | 8.231369 |
| 2012 | 8.081428 |
| 2013 | 8.537331 |
| 2015 | 7.219164 |
| 2016 | 7.791453 |
| 2017 | 7.322790 |
| 2018 | 7.113402 |
| 2019 | 7.880447 |
| 2020 | 8.046531 |
| 2021 | 8.095523 |
| 2022 | 7.621146 |
ggplot(standard_dev_all_431, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 431 average temperatures for water years 2005-2021
sd_mk_431 <- mk.test(standard_dev_all_431$sd_2)
print(sd_mk_431)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -104.0000000 2842.0000000 -0.2561576
sd_sens_431 <- sens.slope(standard_dev_all_431$sd_2)
print(sd_sens_431)
##
## Sen's slope
##
## data: standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.041299369 0.001003281
## sample estimates:
## Sen's slope
## -0.01890302
#using the clean culled df:
#average water year temperature
yearly_wy_aver_431_ad <- snotel_431_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_431_ad <- yearly_wy_aver_431_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_431_ad <- daily_wy_aver_431_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_431_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_431_ad <-daily_wy_aver_431_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_431_ad$date_temp_ad <- signif(daily_wy_aver2_431_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_431_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_431_ad <- daily_wy_aver_431_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431_ad <- standard_dev_431_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_431_ad <- standard_dev_all_431_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_431_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1989 | 7.914872 |
| 1993 | 7.518415 |
| 1995 | 7.173944 |
| 1996 | 7.365597 |
| 1997 | 7.634730 |
| 1998 | 7.970563 |
| 1999 | 6.682639 |
| 2000 | 7.537564 |
| 2001 | 8.039524 |
| 2002 | 8.344015 |
| 2003 | 7.813684 |
| 2004 | 7.737034 |
| 2005 | 7.599212 |
| 2006 | 7.418573 |
| 2007 | 7.938448 |
| 2008 | 8.260644 |
| 2009 | 7.298209 |
| 2010 | 8.419846 |
| 2011 | 8.232506 |
| 2012 | 8.082571 |
| 2013 | 8.538590 |
| 2015 | 7.220737 |
| 2016 | 7.793079 |
| 2017 | 7.324026 |
| 2018 | 7.114853 |
| 2019 | 7.881605 |
| 2020 | 8.047598 |
| 2021 | 8.096826 |
| 2022 | 7.622159 |
ggplot(standard_dev_all_431_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 431 average temperatures for water years 1986-2021
sd_mk_431_ad <- mk.test(standard_dev_all_431_ad$sd_2)
print(sd_mk_431_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 54.0000000 2842.0000000 0.1330049
sd_sens_431_ad <- sens.slope(standard_dev_all_431_ad$sd_2)
print(sd_sens_431_ad)
##
## Sen's slope
##
## data: standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01334615 0.03134200
## sample estimates:
## Sen's slope
## 0.01029449
Morrisey 7/12/2004
snotel_538 <- SNOTEL_san_juan_Area %>%
filter(site_id == "538")
#str(snotel_538) # check the date, usually a character.
snotel_538$Date <- as.Date(snotel_538$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_538_clean <- snotel_538 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_538_clean <- snotel_538_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_538_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_538_clean <- snotel_538_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_min < 25)
ggplot(snotel_538_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_538_cull_count <- snotel_538_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_538_cull_count
## # A tibble: 36 x 2
## # Groups: waterYear [36]
## waterYear n
## <dbl> <int>
## 1 1987 355
## 2 1988 366
## 3 1989 364
## 4 1990 365
## 5 1991 337
## 6 1992 358
## 7 1993 346
## 8 1994 351
## 9 1995 364
## 10 1996 365
## # ... with 26 more rows
# filtering for too few observations in a year
snotel_538_cull_count_days <- snotel_538_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_538_cull_count_days
## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 1991 337
## 2 1993 346
1991 and 1993 need to be culled.
snotel_538_clean_culled <- snotel_538_clean %>%
filter(waterYear != "1991" & waterYear != "1993") # & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_538_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)
dygraph(temp_538_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_538_clean_culled <- snotel_538_clean_culled %>%
filter(temperature_mean < 20)
temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)
dygraph(temp_538_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
Morrisey 7/12/2004
snotel_538_adjusted <- snotel_538_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-07-12", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_538 <- snotel_538_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_538 <- yearly_wy_aver_538 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_538 <- daily_wy_aver_538 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_538$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_538 <-daily_wy_aver_538 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_538$date_temp <- signif(daily_wy_aver2_538$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_538, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_538 <- daily_wy_aver_538 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538 <- standard_dev_538 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_538 <- standard_dev_all_538 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_538 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.006364 |
| 1988 | 8.601864 |
| 1989 | 8.632529 |
| 1990 | 8.392986 |
| 1992 | 7.829013 |
| 1994 | 10.169815 |
| 1995 | 7.658097 |
| 1996 | 7.925011 |
| 1997 | 7.934668 |
| 1998 | 8.455869 |
| 1999 | 7.252626 |
| 2000 | 7.992580 |
| 2001 | 8.428201 |
| 2002 | 8.883209 |
| 2003 | 8.309823 |
| 2004 | 8.066471 |
| 2005 | 7.228205 |
| 2006 | 7.518896 |
| 2007 | 7.890984 |
| 2008 | 8.117130 |
| 2009 | 7.048527 |
| 2010 | 8.062733 |
| 2011 | 7.859864 |
| 2012 | 7.787906 |
| 2013 | 8.410412 |
| 2014 | 7.465408 |
| 2015 | 6.906228 |
| 2016 | 7.611126 |
| 2017 | 7.235848 |
| 2018 | 7.239780 |
| 2019 | 7.718304 |
| 2020 | 7.976082 |
| 2021 | 8.027478 |
| 2022 | 7.664387 |
ggplot(standard_dev_all_538, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 538 average temperatures for water years 2005-2021
sd_mk_538 <- mk.test(standard_dev_all_538$sd_2)
print(sd_mk_538)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -181.0000000 4550.3333333 -0.3226381
sd_sens_538 <- sens.slope(standard_dev_all_538$sd_2)
print(sd_sens_538)
##
## Sen's slope
##
## data: standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.042981556 -0.006896251
## sample estimates:
## Sen's slope
## -0.02384088
#using the clean culled df:
#average water year temperature
yearly_wy_aver_538_ad <- snotel_538_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_538_ad <- yearly_wy_aver_538_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_538_ad <- daily_wy_aver_538_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_538_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_538_ad <-daily_wy_aver_538_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_538_ad$date_temp_ad <- signif(daily_wy_aver2_538_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_538_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_538_ad <- daily_wy_aver_538_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538_ad <- standard_dev_538_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_538_ad <- standard_dev_all_538_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_538_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 7.447584 |
| 1988 | 8.019241 |
| 1989 | 8.069873 |
| 1990 | 7.812557 |
| 1992 | 7.272544 |
| 1994 | 9.549437 |
| 1995 | 7.113641 |
| 1996 | 7.353811 |
| 1997 | 7.396458 |
| 1998 | 7.860736 |
| 1999 | 6.739299 |
| 2000 | 7.408528 |
| 2001 | 7.835100 |
| 2002 | 8.260161 |
| 2003 | 7.706375 |
| 2004 | 7.450714 |
| 2005 | 7.228187 |
| 2006 | 7.519259 |
| 2007 | 7.891254 |
| 2008 | 8.117349 |
| 2009 | 7.048722 |
| 2010 | 8.062844 |
| 2011 | 7.860215 |
| 2012 | 7.788573 |
| 2013 | 8.411054 |
| 2014 | 7.465688 |
| 2015 | 6.906947 |
| 2016 | 7.611742 |
| 2017 | 7.235817 |
| 2018 | 7.240024 |
| 2019 | 7.718445 |
| 2020 | 7.976612 |
| 2021 | 8.027693 |
| 2022 | 7.664991 |
ggplot(standard_dev_all_538_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 538 average temperatures for water years 1986-2021
sd_mk_538_ad <- mk.test(standard_dev_all_538_ad$sd_2)
print(sd_mk_538_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 1.100000e+01 4.550333e+03 1.960784e-02
sd_sens_538_ad <- sens.slope(standard_dev_all_538_ad$sd_2)
print(sd_sens_538_ad)
##
## Sen's slope
##
## data: standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01508416 0.01703536
## sample estimates:
## Sen's slope
## 0.0009105465
Oyler -> Morrisey 6/22/2005
snotel_589 <- SNOTEL_san_juan_Area %>%
filter(site_id == "589")
#str(snotel_589) # check the date, usually a character.
snotel_589$Date <- as.Date(snotel_589$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_589_clean <- snotel_589 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_589_clean <- snotel_589_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_589_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_589_clean <- snotel_589_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_min < 25)
ggplot(snotel_589_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_589_cull_count <- snotel_589_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_589_cull_count
## # A tibble: 37 x 2
## # Groups: waterYear [37]
## waterYear n
## <dbl> <int>
## 1 1986 1
## 2 1987 358
## 3 1988 366
## 4 1989 362
## 5 1990 365
## 6 1991 359
## 7 1992 366
## 8 1993 350
## 9 1994 337
## 10 1995 364
## # ... with 27 more rows
# filtering for too few observations in a year
snotel_589_cull_count_days <- snotel_589_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_589_cull_count_days
## # A tibble: 10 x 2
## # Groups: waterYear [10]
## waterYear n
## <dbl> <int>
## 1 1986 1
## 2 1994 337
## 3 1997 346
## 4 2006 317
## 5 2009 314
## 6 2010 332
## 7 2011 309
## 8 2012 322
## 9 2013 331
## 10 2014 335
1986, 1994, 1997, 2006, 2009, 2010, 2011, 2012, 2013, 2014 need to be culled.
snotel_589_clean_culled <- snotel_589_clean %>%
filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1997" & waterYear != "2006" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_589_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)
dygraph(temp_589_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_589_clean_culled <- snotel_589_clean_culled %>%
filter(temperature_mean < 20)
temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)
dygraph(temp_589_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
Oyler -> Morrisey 6/22/2005
snotel_589_adjusted <- snotel_589_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_589 <- snotel_589_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_589 <- yearly_wy_aver_589 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_589 <- daily_wy_aver_589 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_589$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_589 <-daily_wy_aver_589 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_589$date_temp <- signif(daily_wy_aver2_589$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_589, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_589 <- daily_wy_aver_589 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589 <- standard_dev_589 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_589 <- standard_dev_all_589 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_589 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.224944 |
| 1988 | 8.773281 |
| 1989 | 8.776693 |
| 1990 | 8.580289 |
| 1991 | 8.502502 |
| 1992 | 7.975611 |
| 1993 | 8.199081 |
| 1995 | 7.946318 |
| 1996 | 8.231469 |
| 1998 | 8.706109 |
| 1999 | 7.431979 |
| 2000 | 8.232988 |
| 2001 | 8.771085 |
| 2002 | 9.082003 |
| 2003 | 8.732769 |
| 2004 | 8.380419 |
| 2005 | 8.458900 |
| 2007 | 8.212611 |
| 2008 | 8.408220 |
| 2015 | 7.103666 |
| 2016 | 7.842413 |
| 2017 | 7.527960 |
| 2018 | 7.506884 |
| 2019 | 8.029199 |
| 2020 | 8.271099 |
| 2021 | 8.278625 |
| 2022 | 8.020852 |
ggplot(standard_dev_all_589, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 589 average temperatures for water years 2005-2021
sd_mk_589 <- mk.test(standard_dev_all_589$sd_2)
print(sd_mk_589)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -89.0000000 2301.0000000 -0.2535613
sd_sens_589 <- sens.slope(standard_dev_all_589$sd_2)
print(sd_sens_589)
##
## Sen's slope
##
## data: standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.043902404 0.001230005
## sample estimates:
## Sen's slope
## -0.02183399
#using the clean culled df:
#average water year temperature
yearly_wy_aver_589_ad <- snotel_589_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_589_ad <- yearly_wy_aver_589_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_589_ad <- daily_wy_aver_589_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_589_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_589_ad <-daily_wy_aver_589_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_589_ad$date_temp_ad <- signif(daily_wy_aver2_589_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_589_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_589_ad <- daily_wy_aver_589_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589_ad <- standard_dev_589_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_589_ad <- standard_dev_all_589_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_589_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 7.671688 |
| 1988 | 8.200568 |
| 1989 | 8.222864 |
| 1990 | 7.999469 |
| 1991 | 7.962227 |
| 1992 | 7.420415 |
| 1993 | 7.639311 |
| 1995 | 7.378994 |
| 1996 | 7.641722 |
| 1998 | 8.086353 |
| 1999 | 6.914950 |
| 2000 | 7.633161 |
| 2001 | 8.156564 |
| 2002 | 8.440080 |
| 2003 | 8.103839 |
| 2004 | 7.827141 |
| 2005 | 7.781158 |
| 2007 | 8.212720 |
| 2008 | 8.408841 |
| 2015 | 7.104293 |
| 2016 | 7.842763 |
| 2017 | 7.527954 |
| 2018 | 7.507329 |
| 2019 | 8.029667 |
| 2020 | 8.271310 |
| 2021 | 8.278726 |
| 2022 | 8.021229 |
ggplot(standard_dev_all_589_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 589 average temperatures for water years 1986-2021
sd_mk_589_ad <- mk.test(standard_dev_all_589_ad$sd_2)
print(sd_mk_589_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 3.300000e+01 2.301000e+03 9.401709e-02
sd_sens_589_ad <- sens.slope(standard_dev_all_589_ad$sd_2)
print(sd_sens_589_ad)
##
## Sen's slope
##
## data: standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01153608 0.02498422
## sample estimates:
## Sen's slope
## 0.008250776
Morrisey 6/26/2006
snotel_624 <- SNOTEL_san_juan_Area %>%
filter(site_id == "624")
#str(snotel_624) # check the date, usually a character.
snotel_624$Date <- as.Date(snotel_624$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_624_clean <- snotel_624 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_624_clean <- snotel_624_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_624_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_624_clean <- snotel_624_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)
ggplot(snotel_624_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_624_cull_count <- snotel_624_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_624_cull_count
## # A tibble: 43 x 2
## # Groups: waterYear [43]
## waterYear n
## <dbl> <int>
## 1 1980 59
## 2 1981 354
## 3 1982 1
## 4 1983 312
## 5 1984 360
## 6 1985 362
## 7 1986 364
## 8 1987 360
## 9 1988 340
## 10 1989 363
## # ... with 33 more rows
# filtering for too few observations in a year
snotel_624_cull_count_days <- snotel_624_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_624_cull_count_days
## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1980 59
## 2 1982 1
## 3 1983 312
## 4 1988 340
## 5 1991 338
## 6 1993 330
## 7 1994 330
## 8 2021 338
## 9 2022 348
1980, 1982, 1983, 1988, 1991, 1993, 1994, 2021, 2022 need to be culled.
snotel_624_clean_culled <- snotel_624_clean %>%
filter(waterYear != "1980" & waterYear != "1981" & waterYear != "1982" & waterYear != "1983" & waterYear != "1988" & waterYear != "1991" & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022")# & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_624_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_624_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)
dygraph(temp_624_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_624_clean_culled <- snotel_624_clean_culled %>%
filter(temperature_mean < 20)
temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)
dygraph(temp_624_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_624_adjusted <- snotel_624_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-06-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_624 <- snotel_624_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_624 <- yearly_wy_aver_624 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_624 <- daily_wy_aver_624 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_624$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_624 <-daily_wy_aver_624 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_624$date_temp <- signif(daily_wy_aver2_624$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_624, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_624 <- daily_wy_aver_624 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624 <- standard_dev_624 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_624 <- standard_dev_all_624 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_624 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1984 | 8.239005 |
| 1985 | 8.207104 |
| 1986 | 7.191318 |
| 1987 | 7.871848 |
| 1989 | 8.454185 |
| 1990 | 8.200220 |
| 1992 | 6.438238 |
| 1995 | 7.565327 |
| 1996 | 7.835558 |
| 1997 | 7.996219 |
| 1998 | 8.415139 |
| 1999 | 7.174609 |
| 2000 | 7.839800 |
| 2001 | 8.355597 |
| 2002 | 8.492066 |
| 2003 | 8.150297 |
| 2004 | 7.912597 |
| 2005 | 7.915951 |
| 2006 | 7.936724 |
| 2007 | 7.640492 |
| 2008 | 8.021623 |
| 2009 | 7.171130 |
| 2010 | 8.259271 |
| 2011 | 7.932618 |
| 2012 | 7.791063 |
| 2013 | 8.259320 |
| 2014 | 7.372950 |
| 2015 | 6.806548 |
| 2016 | 7.548199 |
| 2017 | 7.108731 |
| 2018 | 7.036935 |
| 2019 | 7.675916 |
| 2020 | 7.741703 |
ggplot(standard_dev_all_624, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 624 average temperatures for water years 2005-2021
sd_mk_624 <- mk.test(standard_dev_all_624$sd_2)
print(sd_mk_624)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -128.0000000 4165.3333333 -0.2424242
sd_sens_624 <- sens.slope(standard_dev_all_624$sd_2)
print(sd_sens_624)
##
## Sen's slope
##
## data: standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.0333661824 -0.0003479011
## sample estimates:
## Sen's slope
## -0.01662145
#using the clean culled df:
#average water year temperature
yearly_wy_aver_624_ad <- snotel_624_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_624_ad <- yearly_wy_aver_624_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_624_ad <- daily_wy_aver_624_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_624_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_624_ad <-daily_wy_aver_624_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_624_ad$date_temp_ad <- signif(daily_wy_aver2_624_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_624_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_624_ad <- daily_wy_aver_624_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624_ad <- standard_dev_624_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_624_ad <- standard_dev_all_624_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_624_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1984 | 7.738801 |
| 1985 | 7.713059 |
| 1986 | 6.725366 |
| 1987 | 7.374899 |
| 1989 | 7.942209 |
| 1990 | 7.687289 |
| 1992 | 6.080980 |
| 1995 | 7.080976 |
| 1996 | 7.336792 |
| 1997 | 7.519955 |
| 1998 | 7.885653 |
| 1999 | 6.718142 |
| 2000 | 7.329553 |
| 2001 | 7.833709 |
| 2002 | 7.944757 |
| 2003 | 7.617514 |
| 2004 | 7.442397 |
| 2005 | 7.390605 |
| 2006 | 7.304030 |
| 2007 | 7.639661 |
| 2008 | 8.021493 |
| 2009 | 7.170919 |
| 2010 | 8.259153 |
| 2011 | 7.932368 |
| 2012 | 7.790490 |
| 2013 | 8.259010 |
| 2014 | 7.372461 |
| 2015 | 6.806595 |
| 2016 | 7.547948 |
| 2017 | 7.107845 |
| 2018 | 7.036548 |
| 2019 | 7.675791 |
| 2020 | 7.741347 |
ggplot(standard_dev_all_624_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 624 average temperatures for water years 1986-2021
sd_mk_624_ad <- mk.test(standard_dev_all_624_ad$sd_2)
print(sd_mk_624_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 2.400000e+01 4.165333e+03 4.545455e-02
sd_sens_624_ad <- sens.slope(standard_dev_all_624_ad$sd_2)
print(sd_sens_624_ad)
##
## Sen's slope
##
## data: standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01304163 0.02223675
## sample estimates:
## Sen's slope
## 0.003307866
NSCE-Morrisey, Bias-Original 10/2/2003
snotel_632 <- SNOTEL_san_juan_Area %>%
filter(site_id == "632")
#str(snotel_632) # check the date, usually a character.
snotel_632$Date <- as.Date(snotel_632$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_632_clean <- snotel_632 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_632_clean <- snotel_632_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_632_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_632_clean <- snotel_632_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
#filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)
ggplot(snotel_632_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_632_cull_count <- snotel_632_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_632_cull_count
## # A tibble: 37 x 2
## # Groups: waterYear [37]
## waterYear n
## <dbl> <int>
## 1 1986 54
## 2 1987 355
## 3 1988 362
## 4 1989 364
## 5 1990 365
## 6 1991 365
## 7 1992 366
## 8 1993 303
## 9 1994 294
## 10 1995 361
## # ... with 27 more rows
# filtering for too few observations in a year
snotel_632_cull_count_days <- snotel_632_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_632_cull_count_days
## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1986 54
## 2 1993 303
## 3 1994 295
## 4 2005 344
## 5 2006 304
## 6 2013 248
1986, 1993, 1994, 2005, 2006, 2013 need to be culled.
snotel_632_clean_culled <- snotel_632_clean %>%
filter(waterYear != "1986" & waterYear != "1993" & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_632_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_632_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)
dygraph(temp_632_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_632_clean_culled <- snotel_632_clean_culled %>%
filter(temperature_mean < 20)
temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)
dygraph(temp_632_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
NSCE-Morrisey, Bias-Original 10/2/2003
snotel_632_adjusted <- snotel_632_clean_culled %>%
mutate(temp_ad = if_else(Date < "2003-10-02", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_632 <- snotel_632_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_632 <- yearly_wy_aver_632 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_632 <- daily_wy_aver_632 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_632$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_632 <-daily_wy_aver_632 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_632$date_temp <- signif(daily_wy_aver2_632$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_632, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_632 <- daily_wy_aver_632 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632 <- standard_dev_632 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_632 <- standard_dev_all_632 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_632 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 7.613234 |
| 1988 | 8.191410 |
| 1989 | 8.166070 |
| 1990 | 7.646926 |
| 1991 | 7.969691 |
| 1992 | 7.156808 |
| 1995 | 7.489888 |
| 1996 | 7.655841 |
| 1997 | 7.698012 |
| 1998 | 8.157837 |
| 1999 | 6.885991 |
| 2000 | 7.657554 |
| 2001 | 6.788385 |
| 2002 | 8.254532 |
| 2003 | 7.907020 |
| 2004 | 7.308264 |
| 2007 | 7.607147 |
| 2008 | 7.964946 |
| 2009 | 6.957483 |
| 2010 | 7.983696 |
| 2011 | 7.701623 |
| 2012 | 7.490054 |
| 2014 | 7.236529 |
| 2015 | 6.581104 |
| 2016 | 7.409875 |
| 2017 | 7.001441 |
| 2018 | 6.814621 |
| 2019 | 7.582515 |
| 2020 | 7.766470 |
| 2021 | 7.792571 |
| 2022 | 7.357711 |
ggplot(standard_dev_all_632, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 632 average temperatures for water years 2005-2021
sd_mk_632 <- mk.test(standard_dev_all_632$sd_2)
print(sd_mk_632)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -97.0000000 3461.6666667 -0.2086022
sd_sens_632 <- sens.slope(standard_dev_all_632$sd_2)
print(sd_sens_632)
##
## Sen's slope
##
## data: standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.034372524 0.003422909
## sample estimates:
## Sen's slope
## -0.01583493
#using the clean culled df:
#average water year temperature
yearly_wy_aver_632_ad <- snotel_632_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_632_ad <- yearly_wy_aver_632_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_632_ad <- daily_wy_aver_632_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_632_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_632_ad <-daily_wy_aver_632_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_632_ad$date_temp_ad <- signif(daily_wy_aver2_632_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_632_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_632_ad <- daily_wy_aver_632_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632_ad <- standard_dev_632_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_632_ad <- standard_dev_all_632_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_632_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 7.106416 |
| 1988 | 7.678134 |
| 1989 | 7.662133 |
| 1990 | 7.141729 |
| 1991 | 7.492618 |
| 1992 | 6.671512 |
| 1995 | 6.957618 |
| 1996 | 7.110780 |
| 1997 | 7.178235 |
| 1998 | 7.586041 |
| 1999 | 6.397699 |
| 2000 | 7.103150 |
| 2001 | 6.358660 |
| 2002 | 7.672831 |
| 2003 | 7.338249 |
| 2004 | 7.309361 |
| 2007 | 7.606845 |
| 2008 | 7.965266 |
| 2009 | 6.957636 |
| 2010 | 7.984148 |
| 2011 | 7.701741 |
| 2012 | 7.490245 |
| 2014 | 7.236308 |
| 2015 | 6.581364 |
| 2016 | 7.409996 |
| 2017 | 7.000722 |
| 2018 | 6.815246 |
| 2019 | 7.583475 |
| 2020 | 7.766635 |
| 2021 | 7.793272 |
| 2022 | 7.358050 |
ggplot(standard_dev_all_632_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 632 average temperatures for water years 1986-2021
sd_mk_632_ad <- mk.test(standard_dev_all_632_ad$sd_2)
print(sd_mk_632_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 57.0000000 3461.6666667 0.1225806
sd_sens_632_ad <- sens.slope(standard_dev_all_632_ad$sd_2)
print(sd_sens_632_ad)
##
## Sen's slope
##
## data: standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01086011 0.02505935
## sample estimates:
## Sen's slope
## 0.008387789